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Abstract 


A  theory  for  time  correlation  functions  in  liquids  is  developed  based  on  the 
optimized  quadratic  approximation  for  liquid  state  potential  energy  functions. 
The  latter  approximation  leads  to  the  rigorous  mathematical  definition  of  in¬ 
herent  structures  in  liquids  and  their  vibrational  fluctuations)  in  turn  leading 
to  the  concept  of  inherent  normal  modes  in  the  liquid  state.  These  normal 
modes  are  called  “optimized  normal  modes”.  Unlike  normal  modes  based 
on  instantaneous  liquid  state  configurations,  the  optimized  normal  modes  are 


stable,  having  real-valued  frequencies,  and  each  inherent  liquid  state  struc¬ 
ture  has  a  different  set  of  modes  associated  with  it.  By  including  a  single 
phenomenological  decay  function  which  captures  the  average  transition  rate 
between  the  different  sets  of  normal  modes,  velocity  time  correlation  func¬ 
tions  and  dynamical  friction  kernels  for  solute  bonds  can  be  predicted  in 
good  agreement  with  direct  molecular  dynamics  simulation  results. 
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I.  INTRODUCTION 


One  of  the  basic  ingredients  in  condensed  matter  theory  [1,2]  is  the  concept  of  phonons, 
i.e.,  the  small  oscillations  about  a  stable  structure  or  the  energy  minimum,  which  can  be 
related  to  many  equilibrium  and  transport  properties  such  as  heat  capacity,  thermal  con¬ 
ductivity,  thermal  expansion,  light  scattering,  etc.  It  is  clear,  however,  that  liquids  are 
different  from  solids  because  of  their  lack  of  stable  structures,  making  it  formally  difficult 
to  apply  the  well-developed  theory  of  phonons  in  order  to  calculate,  e.g.,  time  correlation 
functions  and  thereby  predict  experimental  observables.  A  clear  challenge,  therefore,  is  to 
derive  from  first  principles  a  set  of  '‘modes”  which  to  some  degree  dominate  the  molecular 
motions  in  liquids,  at  least  for  times  smaller  than  some  phenomenological  relaxation  time  of 
such  modes.  [3,4]  The  purpose  of  this  paper  is  to  provide  a  formal  prescription  for  defining 
such  modes  which  we  call  “optimized  normal  modes”  (ONM).  A  related  and  challenging 
problem  is  the  microscopic  origin  of  the  relaxation  behavior  of  these  modes,  but  this  issue 
will  be  left  to  future  research. 

It  should  be  noted  that  the  theory  described  herein  builds  on  the  earlier  work  of  Stillinger 
and  Weber  [5-7]  and  of  Zwanzig.  [8]  The  former  authors  proposed  the  inherent  structure 
picture  of  liquids  in  which  such  structures  are  determined  by  a  steepest-decent  quench  on 
the  liquid  state  potential  hypersurface.  The  many-body  phase  space  is  thus  divided  into 
subspaces  corresponding  to  many  different  local  minima.  The  distribution  of  the  inherent 
structure  local  minima  depends  on  the  interaction  potential,  temperature,  and  density.  The 
Stillinger  and  Weber  stable  states,  [5-7]  being  local  potential  minima,  are  free  of  imaginary 
frequencies  and  thus  ideal  candidates  for  an  effective  harmonic  approximation.  Unlike  in 
solids,  however,  the  inherent  structure  is  a  metastable  state  so  there  must  be  an  overall 
decay  behavior  associated  with  the  transitions  between  the  metastable  states.  Though  rich 
in  physical  insight,  the  work  by  Stillinger  and  Weber  did  not  provide  a  variational  procedure 
for  defining  the  inherent  structures  and  their  associated  vibrational  modes  for  a  given  set  of 
thermodynamical  state  variables  (e.g.,  temperature). 
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In  his  dynamical  view  of  liquids,  Zwanzig  [8]  suggested  that  after  a  period  of  vibrational 
motion  about  an  inherent  structure  minimum,  the  liquid  will  jump  through  a  saddle  point 
to  another  local  minimum  and  its  associated  vibrations.  Eventually,  the  liquid  will  explore 
all  the  phase  space  available  to  it.  The  transition  process  is  characterized  by  an  average 
lifetime  “r”  and  thus  an  exponential  decay  factor  is  imposed  on  the  harmonic  motion. 
The  identification  of  the  inherent  structures,  the  vibrational  motion  about  them,  and  the 
inter-minima  transitions  provides  a  plausible  picture  of  the  underlying  dynamical  behavior 
of  a  liquid.  Based  on  this  picture,  Zwanzig  [8]  derived  an  expression  for  the  self-diffusion 
constant  and  a  relation  between  self-diffusion  and  viscosity  which  is  consistent  with  the 
Stokes-Einstein  law.  In  this  theory,  Zwanzig  invoked  the  idea  of  inherent  structure  modes, 
but  he  did  not  explicitly  define  those  modes  starting  with  the  microscopic  potential  (i.e.,  he 
used  a  Debye-like  approximation).  In  the  present  paper,  a  mathematical  procedure  is  used  to 
specify  the  inherent  structure  modes  (i.e.,  the  optimized  normal  modes)  which  provides  the 
missing  element  in  Zwanzig’s  picture.  (See  also  Refs.  [9]  and  [10]  for  another  such  approach.) 

In  our  previous  study  on  the  formulation  of  statistical  mechanics  based  on  an  effective 
quadratic  potential,  [11]  the  exact  cumulant  expansion  of  the  partition  function  was  shown 
to  have  a  one-to-one  correspondence  with  a  diagrammatic  representation.  It  was  also  shown 
that  diagrammatic  classifications  and  topological  reductions  result  in  the  renormalization  of 
the  three  diagrammatic  elements  and  thus  lead  to  a  set  of  self-consistent  effective  quadratic 
equations  at  different  levels  of  approximation.  The  theory  is  applicable  to  both  classical  and 
quantum  systems,  and  can  be  shown  in  the  extremely  low  temperature  limit  to  be  equivalent 
to  ground  state  calculations  in  a  harmonic  oscillator  basis  set.  Among  the  central  results 
of  the  formalism  [11]  is  the  optimized  quadratic  approximation  (OQA)  for  the  partition 
function  which  is  of  special  importance  because  of  its  applicability  to  realistic  many-body 
systems.  The  lowest-order  OQA  equations  can  be  derived  from  the  diagrammatic  approach 
or,  equivalently,  from  the  Gibbs-Bogoliubov  variational  principle.  It  can  be  shown  that  the 
optimized  quadratic  reference  potential  represents  the  best  fit  to  an  anharmonic  potential 
and  any  further  corrections  are  beyond  the  effective  potential  description. 
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Similar  to  the  inherent  structure  picture,  there  exist  multiple  solutions  to  the  OQA 
equations  at  different  positions,  each  corresponding  to  a  local  potential  minimum  with  its 
associated  effective  normal  mode  frequencies.  Therefore,  the  OQA  solutions  separate  the 
many-body  potential  hypersurface  into  different  regions,  each  having  a  definite  thermal 
partition,  and  all  physical  quantities  can  be  expressed  as  a  superposition  of  the  OQA  ex¬ 
pectation  values  weighted  by  the  thermal  partitions.  The  nature  of  the  OQA  solution  on  a 
many-body  hypersurface  reveals  the  distinction  between  solids,  liquids,  and  glasses,  as  well 
as  the  transitions  between  those  phases.  When  applied  to  liquids,  this  theory  defines  the 
inherent  structures  within  a  rigorous  theoretical  framework  and,  furthermore,  introduces 
the  optimized  normal  modes  (ONM)  of  oscillation  in  a  well-defined  fashion.  Thus,  to  a 
reasonable  degree  of  mathematical  rigor,  there  exist  solid  state-like  concepts  in  liquids  such 
as  equilibrium  structures  and  phonon  excitations,  though  they  are  of  a  metastable  nature. 

The  notion  of  inherent  structures  and  their  optimized  normal  modes  immediately  induces 
one  to  extend  the  OQA  theory  to  characterize  the  dynamics  of  liquid  state  systems,  partic¬ 
ularly  time  correlation  functions.  At  the  level  of  the  OQA,  the  approach  expresses  the  time 
correlation  function  as  a  thermal  partition-weighted  superposition  of  optimized  harmonic 
oscillator  time  correlation  functions.  As  suggested  by  Zwanzig,  [8]  such  a  linear  description 
is  adequate  for  a  time  interval  shorter  than  some  relaxation  time,  beyond  which  the  effective 
harmonic  motion  for  an  inherent  structure  decays  into  the  effective  harmonic  motion  for  a 
different  inherent  structure.  Motivated  by  this  physical  intuition,  a  decay  factor  is  incorpo¬ 
rated  into  the  expression  for  the  time  correlation  function  predicted  by  the  optimized  normal 
modes,  thus  introducing  a  broadening  of  the  spectrum  which  defines  damped  normal  modes 
(DNM).  The  decay  function  describes  the  average  long  time  decay  of  correlations  due  to  the 
transitions  between  the  normal  modes  of  different  inherent  structures.  It  will  be  shown  that 
the  combination  of  the  DNM  picture  with  the  self-consistent  OQA  proves  to  be  a  fruitful 
theoretical  framework  for  predicting  liquid  state  time  correlation  functions. 

Another  outgrowth  of  the  present  theory  applies  to  a  wide  variety  of  processes  involv¬ 
ing  intramolecular  motions  in  liquids  which  can  be  modeled  by  the  generalized  Langevin 
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equation  (GLE).  [3,12,13]  In  the  GLE  approach,  the  dynamical  solvent  effect  on,  e.g.,  a 
molecular  bond  or  some  more  complicated  solute  coordinate  is  characterized  by  a  dynamic 
friction  kernel  which  can  be  predicted  theoretically  only  in  the  simple  case  of  translational 
and  rotational  motions.  [14,15]  Though  the  validity  of  the  GLE  has  not  been  proven  in  gen¬ 
eral,  it  can  be  rigorously  derived  for  a  Gaussian  bath  which  consists  of  harmonic  oscillators 
linearly  coupled  to  the  solute.  [4]  It  will  be  shown  that  the  OQA  and  DNM  theory  can 
be  extended  to  identify  the  bath  normal-modes  in  a  many-body  system  and  their  coupling 
coefficients  to  a  solute,  thus  providing  a  theory  for  the  friction  kernel  in  the  GLE  picture. 
This  generalization  of  the  theory  has  important  implications  for  the  study  of  friction  on 
solute  bond  vibrations  as  well  as  activated  barrier  crossing.  In  accord  with  the  Zwanzig 
picture,  the  GLE  bath  from  the  DNM  theory  experiences  an  exponential  decay  because  of 
the  transitions  between  different  inherent  structure  OQA  solutions.  By  assuming  that  the 
exponential  decay  parameter  can  be  determined  from  the  self-diffusion  behavior  of  the  pure 
solvent,  one  can  predict  the  dynamical  friction  kernel  for  a  solute  bond  in  good  agreement 
with  the  exact  result  from  MD  simulations.  [16]  We  believe  this  result  to  be  significant 
because,  in  principle,  it  provides  a  microscopic  theory  for  the  GLE  friction  kernel  in  liquids. 

Before  proceeding  to  the  next  section,  it  is  important  to  contrast  the  present  theory  with 
the  concept  of  “instantaneous  normal  modes”  (INM).  [9,10,17-25]  In  the  latter  theory,  the 
liquid  state  potential  is  Taylor  expanded  at  different  instantaneous  configurations  through 
quadratic  order.  A  set  of  normal  modes  is  then  obtained  by  diagonalizing  the  force  con¬ 
stant  matrix  and  the  short  time  dynamics  resulting  from  that  liquid  configuration  can  be 
predicted.  This  effective  harmonic  motion  is  suggested  to  persist  up  to  some  characteristic 
relaxation  time,  at  which  point  it  is  transformed  into  motion  characteristic  of  another  set 
of  instantaneous  normal  modes.  [9,10,17-25]  The  overall  short  time  dynamics  of  the  liquid 
is  thereby  determined  by  a  superposition  of  the  harmonic  motions  of  all  possible  configu¬ 
rations.  The  liquid  state  “phonon  spectrum”  is  taken  to  be  the  ensemble  average  of  the 
instantaneous  normal  modes  of  the  liquid  configurations.  [17]  Instantaneous  normal  modes 
have  been  used  to  study,  e.g,  the  short-time  dynamics  of  coupled  translational  and  rotational 
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motions  in  molecular  fluids.  [19]  The  predictions  of  the  short-time  harmonic  motion  were 
compared  with  exact  molecular  dynamics  (MD)  simulation  results  and  found  to  agree  only 
for  short  times.  As  a  result  of  the  anharmonicity  in  the  liquid,  the  difficulty  in  describing 
such  correlation  functions  with  the  INM  theory  arises  due  to  the  presence  of  the  unstable 
modes  which  diverge  exponentially  with  time.  From  this  point  of  view,  it  is  mainly  the  un¬ 
stable  modes  which  destroy  the  linear  motion  of  the  liquid.  Since  the  imaginary  frequencies 
presumably  become  operational  after  the  characteristic  relaxation  time,  it  can  be  argued 
that  the  unstable  modes  should  be  removed  from  the  INM  correlation  function.  [19,20]  Rea¬ 
sonable  agreement  with  exact  MD  correlation  functions  can  be  obtained  when  this  “stable 
mode”  approximation  is  implemented.  [19,20,26,27]  The  underlying  ONM  spectrum  in  the 
DNM  theory  does  not  contain  globally  unstable  modes,  and  it  involves  a  rather  different  set 
of  assumptions  and  approximations  in  its  formulation.  The  predictions  of  the  stable  mode 
INM  and  DNM  theories  will  be  compared  numerically  for  several  examples  in  Sec.  IV. 

The  sections  of  this  paper  are  organized  are  follows:  In  Sec.  II,  the  DNM  analysis  of 
inherent  structures  is  presented  within  the  rigorous  OQA  framework.  Then,  in  Sec.  Ill  the 
GLE  picture  is  formulated  in  terms  of  constrained  OQA  equations  and  the  DNM  theory. 
Some  illustrative  numerical  examples  are  presented  in  Sec.  IV  and  concluding  remarks  are 
given  in  Sec.  V. 

II.  DAMPED  NORMAL  MODE  MODE  THEORY  OF  LIQUIDS 

In  an  effort  to  bring  conceptual  order  to  the  disordered  liquid  state,  Stillinger  and  Weber 
[5-7]  advanced  the  idea  of  separating  the  statistical  mechanical  description  of  liquids  into  two 
distinct  parts,  namely,  the  mechanically  stable  packing  part  and  the  thermally  fluctuating 
part.  Their  key  idea  was  a  configurational  mapping  where  arbitrary  sets  of  molecular  posi¬ 
tions  are  referred  to  potential  minima  which  are  the  inhcvcnt  stvuctuvcs  underlying  the  liquid 
phase.  This  mapping  is  generated  by  a  quench  procedure  which  follows  the  steepest-descent 
paths  on  the  hypersurface  of  the  many-body  system.  Much  attention  has  been  focused  on 
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th.6  identification  and  characterization  of  these  mechanically  stable  packings. 

The  formalism  of  Stillinger  and  Weber  [5]  begins  with  the  canonical  partition  function 
for  N  structureless  classical  particles,  given  by 

Zn  =  j  dqex'p[-pV{q)]  ,  (2-1) 

where  N  is  the  particle  number  and  the  usual  normalization  factor  is  omitted  here  for 
simplicity.  The  configurational  integral  is  next  broken  into  the  separate  contributions  from 
each  quench  region,  i.e., 

Zn  =  [  dqexp[-/?V(q)]  ,  (2-2) 

where  R{1)  defines  the  segment  on  the  potential  hypersurface  which  can  be  uniquely  mapped 
to  the  inherent  structure  designated  by  the  index  1.  Within  the  region  R{1),  any  set  of 
coordinates  can  be  traced  to  the  quenched  inherent  structure,  giving 

V{q)  =  V(q/)  +  AiV^(q)  ,  (2-3) 

where  V{qi)  is  the  potential  local  minimum,  satisfying 

Vy(qi)  =  0  ,  (2-4) 

and  V  denotes  the  spatial  derivative.  Consequently,  the  partition  function  can  be  rewritten 
as 

Zn  =  IIexp(-i3y(qi))  /  dqexp[-/?Azy(q)]  ,  (2-5) 

I  JR{1) 

where  the  integration  accounts  for  the  thermal  fiuctuations  around  the  stable  packing  struc¬ 
ture. 

While  the  quench  procedure  may  reveal  the  hidden  structures  of  the  liquid  phase,  it  may 
not  be  particularly  successful  in  recovering  the  equilibrium  properties  of  liquids  and  even 
less  successful  in  predicting  their  dynamical  properties.  While  the  thermal  fiuctuations  in 
the  inherent  structure  potential  wells  are  suggestive  of  a  linear  harmonic  motion,  it  turns  out 
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that  thermally  broadening  the  quenched  structure  by  using  an  Einstein  or  Debye  vibrational 
approximation  fails  to  reconstruct  important  features  such  as  pair  correlation  functions. 
There  are  at  least  two  reasons  for  this,  one  being  the  significant  anharmonicities  of  the 
liquid  state  interaction  potential  and  the  other  being  the  geometric  disorder  of  the  inherent 
structures.  Thereby,  an  actual  set  of  effective  harmonic  modes  will  bear  little  resemblance 
to  the  phonon  spectrum  of  solids  as  described  by  the  Einstein  or  Debye  models.  Evidently, 
a  systematic  theory  is  required  to  formaUze  the  concept  of  inherent  structures,  to  establish 
the  relationship  between  liquid  dynamics  and  collective  harmonic  motions,  and  to  allow  for 
higher-order  corrections.  This  goal  can  be  accomplished  within  the  framework  of  the  OQA 
theory  developed  for  general  potentials  in  Ref.  [11]- 

To  start,  the  OQA  of  Ref.  [11]  for  an  iV-body  coordinate  space  is  briefly  reviewed.  Here, 
vectors  and  matrices  are  denoted  by  bold  fonts,  and  optimized  frequencies  and  positions  are 
denoted  by  bars.  The  basic  OQA  equations  can  be  written  as 

(Vy(q  +  q))c  =  0  (2.6) 

(V:  Vl/(q-bq))c  =  K  (2.7) 

where  K  is  the  optimized  effective  force  constant  matrix  and  V  is  the  partial  derivative 
vector  Vi  =  d/dqi-  The  notation  {•••)c  here  denotes  a  multidimensional  Gaussian  average 
centered  at  q,  i.e., 

(y(q  +  q))c  =  !  I  ■—  /dql7(q  +  q)  exp[-q-C~^-q/2]  ,  (2.8) 

y'det[27rC] 

where  the  Gaussian  width  factor  matrix  C  is  formally  expressed  in  the  classical  limit  as 

C  =  (/3K)-^  .  (2.9) 

In  terms  of  the  eigensolutions,  a  unitary  matrix  U  can  be  found  to  diagonalize  the  mass- 
scaled  force  constant  matrix  K,  giving  the  effective  normal  modes,  i.e., 

U^KU=[I-(5^]  ,  (2.10) 
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where  {a;}  is  the  set  of  the  eigenfrequencies,  I  is  the  3iV-dimensional  identity  matrix,  and 
[I  •  denotes  a  diagonal  matrix  with  the  i***  diagonal  element  given  by  uf.  The  Gaussian 
width  factor  matrix  in  Eq.  (2.9)  can  also  be  determined  from  the  relation 

C  =  u[l-a]U^  ,  (2.11) 

where  U  =  m  is  the  diagonal  mass  matrix,  and  the  individual  elements  of  the 

normal  mode  thermal  width  vector  are  given  in  the  classical  limit  by 

di  =  l/PCj-f  .  (2.12) 

Thus,  the  set  of  optimized  frequencies  {a;}  and  average  positions  {q}  are  variationally  ob¬ 
tained  as  the  self-consistent  solution  to  the  transcendental  matrix  equations  Eqs.  (2.6)  and 
(2.7)  in  AT’-dimensional  space.  The  quantum  generalization  of  the  OQA  equations  is  given 
in  the  Appendix. 

As  it  stands,  there  are  many  possible  solutions  to  the  self-consistent  OQA  equations, 
each  being  defined  in  a  local  potential  well  of  the  many-body  system.  Under  the  condi¬ 
tion  that  different  wells  are  separated,  i.e.,  the  barrier  between  any  two  neighboring  wells 
is  significantly  higher  than  the  average  thermal  energy,  one  can  assume  a  linear  superpo¬ 
sition  of  all  the  metastable  solutions.  In  this  spirit,  the  partition  function  in  Eq.  (2.1)  is 
intrinsically  separated  into  different  integration  regions  and  can  be  written  in  the  quadratic 
approximation  as  [11] 

Zn  ^  5Z-^i.re/exp[-/?(AU(q/ -t-q))c,]  (2-13) 

I 

where  subscript  I  denotes  the  set  of  OQA  solutions  to  Eqs.  (2.6)  and  (2.7).  In  this 
sense,  the  differences  between  liquids  and  solids  can  be  attributed  to  the  nature  of  the 
OQA  solutions  for  the  many-body  configurations.  Indeed,  this  important  concept  makes 
it  possible  to  rigorously  represent  the  inherent  structures  for  liquids,  which  was  previously 
proposed  and  pursued  from  the  perspective  of  quenched  potential  minima,  and  to  introduce 
an  optimized  normal  mode  spectrum  which  will  be  an  analogue  to  the  phonon  spectrum  of 
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solids.  [1,2]  To  be  more  explicit,  one  can  identify  q,  i.e.,  the  solution  to  Eq.  (2.6),  as  the 
inherent  structures,  and  the  corresponding  eigenvectors  and  eigenvalues,  i.e.,  the  solution  to 
Eq.  (2.7),  as  the  optimized  normal  modes  (ONM)  for  the  inherent  structure.  This  definition 
differs  significantly  from  that  of  Stillinger  and  Weber’s  quenched  minima  for  it  incorporates 
the  packing  structures  and  thermal  fluctuations  into  a  unified  theory.  The  equilibrium  state 
of  the  inherent  structure  defined  in  the  OQA  will  shift  from  the  mechanical  equilibrium  state 
because  of  the  thermal  motion,  while  the  distribution  of  the  optimized  normal  modes  will 
display  very  different  features  from  the  Einstein  or  Debye  model  which  are  only  meaningful 
for  well-defined  solid  lattices.  Moreover,  the  formulation  here  is  applicable  to  both  classical 
and  quantum  Boltzmann  liquids  (cf.  the  Appendix).  Along  these  lines,  we  note  that  the 
mechanical  equilibrium  state  in  Stillinger  and  Weber’s  theory  lacks  a  plausible  interpretation 
in  quantum  mechanics  because  inherent  quantum  fluctuations  must  introduce  uncertainties 
in  the  particle  positions. 

In  general,  all  expectation  values  of  physical  variables  can  be  expressed  as  the  sum  of 
the  distinct  OQA  solutions  weighted  by  the  partitions  of  the  metastable  wells.  To  take  into 
account  the  weight  of  each  solution  correctly,  the  mapping  method,  introduced  as  the  main 
ingredient  in  Stillinger  and  Weber’s  approach,  proves  to  be  helpful.  In  this  approach,  the 
liquid  hypersurface  is  divided  into  regions  R{1),  each  of  which  can  be  mapped  uniquely  to  an 
OQA  inherent  structure.  Correspondingly,  certain  instantaneous  liquid  configurations  are 
traced  to  the  same  OQA  solution,  while  other  OQA  inherent  structures  must  correspond  to 
different  subsurfaces  on  the  liquid  potential  hypersurface.  In  practice,  a  mapping  procedure 
can  be  devised  as  follows: 

(a)  Randomly  select  an  instantaneous  liquid  configuration  from  the  canonical  distribution; 

(b)  Quench  the  liquid  configuration  to  its  potential  minimum; 

(c)  Solve  the  self-consistent  optimized  quadratic  approximation  near  the  mechanical  stable 
structure; 

(d)  Collect  data  for  the  optimized  inherent  structure; 

(e)  Repeat  steps  (a) -(d)  for  many  independent  liquid  configurations. 
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Because  of  the  one-to-one  correspondence  of  the  mapping,  the  thermal  partition  of  the 

optimized  quadratic  solutions  is  accurately  incorporated  into  the  scheme.  The  quench  in 

the  procedure  can  be  achieved  by  the  steepest-descent  or  the  conjugated  gradient  algorithms. 

Many  physical  properties  can  be  investigated  within  the  frame  of  the  optimized  quadratic 

theory.  In  particular,  the  transition  between  the  solid  and  the  liquid  states  can  be  viewed 

as  the  change  from  a  global  inherent  structure  to  many  possible  metastable  structures.  This 

paper,  however,  is  devoted  to  the  study  of  liquid  state  dynciTTiics  based  on  this  model. 

To  begin  the  dynamical  analysis,  one  can  define  an  ONM  spectrum  such  that 

1  ZN 

Donm{'^)  —  2^  (^0)1)90  (2-14) 

where  uit(Qo)  3-^6  the  set  of  eigensolutions  to  Eq.  (2.10)  for  the  region  R{1)  mapped  from 
an  instantaneous  liquid  configuration  qo-  The  average  in  Eq.  (2.14)  is  firom  the  canonical 
distribution  of  liquid  configurations.  The  optimized  normal  modes  variationally  capture 
the  characteristic  stable  mode  thermal  excitations,  [11]  at  least  for  the  time  period  when 
the  system  remains  in  the  metastable  potential  well.  In  contrast,  the  INM  description 
[9,10,17-25]  is  instantaneous  by  nature  and,  through  the  continuity  of  the  fluid  motion, 
renders  instabilities  to  some  of  the  modes.  The  ONM  spectrum  is  essentially  the  liquid  state 
analogue  of  the  self-consistent  phonon  spectrum  of  anharmonic  solids.  [28—30] 

On  the  other  hand,  the  metastability  of  the  optimized  normal  modes  tends  to  destroy 
the  coherence  of  their  vibrations.  The  liquid  motion  can  be  viewed  as  transitions  from 
one  optimized  inherent  structure  to  another,  an  interplay  of  barrier  crossings  and  thermal 
vibrations.  The  short  lifetime  of  the  ONM’s  will  broaden  the  overall  spectrum  as  in  the  case 
of  a  damped  oscillator.  This  dynamical  picture  is  included  by  introducing  a  decay  factor 
into  the  time  correlation  functions.  For  example,  the  velocity  time  correlation  function  for 
a  simple  atomic  fluid  may  be  written  as 

CDNM{t)  =  :^jdujDoNM{^)cos{ujt)f{u},t)  ,  (2.15) 

where  the  subscript  DNM  stand  for  the  “damped  normal  mode”  approximation  and 

is  a  decay  function  which  may,  in  the  most  general  case,  depend  on  the  frequency.  A 


simplifying  assumption  horo  is  to  adopt  a  simplo  monotonic  docay  function  which  ignores 
the  frequency  dependence  such  that  the  DNM  spectrum  now  reads 

Ddnm{(^)  =■  j  du' F{u  —  u')Donm{‘^')  i  (2-16) 

where  F{lj)  is  the  Fourier  transform  of  the  universal  decay  function  which  broadens  the 
ONM  spectrum. 

Since  a  procedure  to  determine  the  functional  form  of  the  decay  function  from  first  prin¬ 
ciples  has  not  yet  been  developed,  it  can  simply  be  assumed  to  be  an  exponential  function, 
i.e., 

f{t)  =  exp(-Altl)  (2.17) 


with  the  corresponding  spectral  convolution  function  F{(jj)  given  by 


F{u)  = 


1  A 

TT  A^  -t- 


(2.18) 


which  takes  the  form  of  a  Lorentzian  broadening  function.  Since  f{t)  is  not  an  analytic 
function  at  t  =  0,  the  decay  constant  A  cannot  be  determined  from  the  short  time  behavior  of 
a  time  correlation  function.  Prom  the  microscopic  point  of  view,  A  can  be  understood  as  the 
average  escape  rate  from  the  OQA  inherent  structures.  Therefore,  one  might  estimate  this 
constant  from  transition  state  theory  (TST)  [13,31]  provided  the  typical  values  of  the  ONM 
frequencies  and  the  barrier  heights  are  available.  Alternatively,  we  note  that  Eq.  (2.16)  with 
the  Lorentzian  broadening  factor  Eq.  (2.18)  is  exactly  the  expression  introduced  by  Zwanzig 
in  his  analysis  of  the  self-diffusion  constant.  [8]  In  the  context  of  the  present  DNM  theory, 
one  can  therefore  adjust  the  decay  constant  A  to  yield  the  correct  experimental  diffusion 
constant,  i.e.. 


D  =  [du  —  r-Y— — ^Donm{^)  )  (2.19) 

J  TT  +UJ^ 

which  is  the  zero-frequency  component  of  the  Fourier  transform  of  the  velocity  autocorrela¬ 
tion  function.  With  the  quantities  Donm{^)  ^  hand,  other  liquid  state  correlation 
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functions  can  then  be  predicted  (cf.  the  next  section).  A  similar  approach  has  been  proposed 
in  Ref.  [10]  in  order  to  introduce  damping  into  the  stable  mode  INM  theory. 

Before  proceeding  to  the  applications  of  the  theory,  we  note  that  an  alternative  choice 
of  f{i)  is  a  Gaussian  decay  function,  i.e., 

f{t)  =  exp(-/ct^/2)  ,  (2.20) 

where  k  is  an  undetermined  constant,  and  the  spectral  convolution  function  is  given  in  this 
case  by 

F{uj)  =  yjexp(-a;^/2K)  (2.21) 

which  broadens  a  single  frequency  into  a  Gaussian  distribution  with  a  width  factor  k.  The 
Gaussian  choice  of  damping  function  f{t)  allows  one  to  impose  the  condition 

(jg  =  J  diij  Donm{!F)oJ^  +  K  (2.22) 

where  is  the  “Einstein  frequency”  calculated  from  equilibrium  properties  via  the  second- 
order  moment  expansion.  This  approach  uniquely  specifies  the  Gaussian  width  factor  k. 
The  optimal  choice  of  the  damping  function  f{t)  will  depend  on  the  problem  at  hand  and 
the  timescale  of  the  behavior  under  examination.  The  two  choices  described  here  are,  of 
course,  qualitatively  different  since  one  (the  Gaussian)  is  based  on  short  time  (inertial) 
behavior,  while  the  other  (the  exponential)  is  based  on  longer  time  (diffusive)  behavior.  In 
the  examples  studied  in  Sec.  IVB,  an  exponential  damping  function  was  found  to  be  superior, 
but  this  is  not  necessarily  always  true. 

III.  DYNAMICAL  FRICTION  ON  SOLUTE  BONDS 

The  generalized  Langevin  equation  (GLE)  has  been  used  to  understand  a  wide  range 
of  problems  involving  molecular  motion  in  liquids  such  as  activated  barrier  crossing  and 
vibrational  relaxation.  [3,4,12,13]  The  GLE  for  a  molecular  bond  can  be  expressed  as 
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(3.1) 


mq(t)  +  mufqit)  +  /  dt'ri{t  -  i)q{i)  =  F(f] 

J  0 

where  q  is  the  displacement  in  the  bond  length,  F{t)  is  the  random  force  along  the  bond,  and 
u  can  be  determined  by  the  mean  square  displacement  =  [/?m(?^)].  Using  projection 
operator  techniques,  one  can  explicitly  derive  the  expressions  for  the  dynamic  friction  rj{t), 
the  random  force  F{t),  and  the  second  fluctuation  dissipation  theorem  which  relates  the 
two.  However  the  formal  definitions  provide  little  help  in  evaluating  these  quantities.  It  is 
thus  necessary  to  obtain  the  dynamical  friction  kernel  by  some  other  means.  An  often  used 
approximation  is  to  set  the  dynamical  friction  kernel  equal  to  the  autocorrelation  function 
of  the  fluctuating  force  exerted  on  the  rigid  bond  by  the  bath  degrees  of  freedom.  The 
rigid  bond  approximation  has  been  shown  to  be  the  high  frequency  limit  of  true  dynamical 
friction  coefficient.  [16] 

If  the  bond  motion  can  be  characterized  by  a  high  frequency  oscillation,  the  dynamic 
friction  kernel  is  equivalent  to  that  evaluated  for  a  rigid  bond  fixed  at  the  average  position  of 
the  bond  coordinate.  [16]  Then,  the  second  dissipation  theorem  yields  a  simple  prescription 
for  the  friction,  i.e., 

r]{t)  =  l3{6F{t)5F{0)),  (3.2) 

where  the  random  force  fluctuation  6F{t)  is  evaluated  with  the  bond  frozen  at  its  equilibrium 
length.  The  explicit  relation  between  the  force  fluctuations  and  the  dynamic  friction  cannot 
be  derived  in  general  except  if  the  nonlinear  bond  coordinate  is  bilinearly  coupled  to  a 
harmonic  bath,  i.e., 

V  =  Veg{q) 

where  Veq{q)  is  the  potential  of  mean  force  along  q,  Xi  is  the  i***  Gaussian  bath  normal  mode, 
and  c,  is  the  coupling  strength.  It  was  shown  by  Zwanzig  [4]  that  the  elimination  of  the  bath 
modes  from  the  equations  of  motion  for  the  above  potential  yields  the  GLE.  The  dynamical 
friction  coefficient  is  then  identified  as 


iV  _2  9  roo 

ri(t)  =  cos{u}it)  =  -  /  du; - cos{ujt)  ,  (3.4) 

^  '  ..  LO-  TT  JO  a; 

1=1 

where  J(a;)  is  the  spectral  density,  defined  in  the  discrete  limit  by 

■  (3-5) 

^  i=l 

The  random  force  can  be  explicitly  expressed  in  terms  of  the  initial  conditions  of  the  bath 
variables.  Therefore,  under  the  assumption  that  the  initial  bath  distribution  in  phase  space 
is  in  thermal  equilibrium  in  the  presence  of  the  system,  one  can  readily  show  that 

,(t)  =  /?(F(()F(0))  ,  (3.6) 

where  the  equilibrium  condition  (F)  =  0  is  implied. 

The  introduction  of  the  spectral  density  J{uS)  makes  it  possible  to  pass  from  a  discrete 
set  of  modes  to  a  continuum  spectrum,  and  hence  to  represent  an  arbitrary  time  dependent 
friction  T]{t^.  The  relation  in  Eq.  (3.6)  holds  for  a  harmonic  bath  regardless  of  the  anharmonic 
bond  potential  or  the  bond  length.  It  is  for  this  reason  that  the  Gaussian  bath  is  an 
attractive  analytical  model  to  study  the  solvent  frictional  effects  on  vibrational  relaxation 
and  activated  reaction  dynamics.  In  a  real  system,  the  asymptotic  limit  of  the  friction 

mentioned  previously  implies  that  the  frequency  of  the  oscillating  bond  must  be  much  larger 

than  the  peak  frequency  of  the  solvent  spectral  density. 

While  the  GLE  is  an  appealing  picture,  questions  remain  whether  the  harmonic  bath  is 
suitable  for  describing  realistic  systems  and,  if  so,  how  the  spectral  density  can  be  calculated 
from  first  principles.  The  rigor  of  such  a  derivation  relies  on  the  quadratic  nature  of  the 
bath  and  the  linearity  of  the  couplings.  This  situation  may  be  best  realized  in  the  solid  state 
where  the  bath  modes  can  be  well  understood  as  the  phonon  excitations.  In  contrast,  there 
is  not  such  a  global  quadratic  bath  for  liquids,  a  bath  defined  as  being  independent  of  the 
temperature  and  density.  In  liquids,  however,  one  can  appeal  to  the  concept  of  an  effective 
Gaussian  bath— precisely  the  target  of  this  paper! 

In  the  previous  section,  it  was  proposed  that  the  configuration  space  of  liquids  can  be 
partitioned  into  different  optimized  metastable  potential  subspaces  so  that  the  short  time 
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liquid  motion  is  described  as  effective  harmonic  thermal  excitations  in  the  inherent  structure 
wells.  The  optimized  structure  and  the  effective  thermal  fluctuations  can  then  be  found  in  a 
self-consistent  fashion  by  virtue  of  the  general  OQA  theory.  Furthermore,  the  dynamics  of 
liquids  can  be  separated  into  the  effective  harmonic  oscillations  in  the  potential  wells  and  the 
transitions  between  the  different  wells.  In  this  simplified  picture,  liquids  can  be  described 
as  a  set  of  damped  harmonic  oscillators.  Following  the  same  line  of  reasoning,  one  can 
picture  a  solute  bond  in  a  solvent  as  being  coupled  to  a  damped  harmonic  bath  consisting  of 
exponentially  decaying  normal  modes.  In  light  of  the  present  theoretical  developments,  the 
harmonic  bath  can  be  identified  as  a  set  of  optimized  normal  modes  in  the  OQA  theory.  The 
normal  modes  thus  defined  will  depend  on  the  particular  liquid  configuration  from  which 
the  optimized  configuration  is  mapped.  Therefore,  the  actual  modes  and  couplings  not  only 
depend  on  the  nature  of  the  solvent,  but  also  on  the  mass,  flexibility,  and  length  of  the  bond. 

In  order  to  formalize  the  DNM  picture  of  dynamical  friction,  a  modification  of  the  OQA 
equations  is  necessary:  the  ONM  solutions  will  apply  to  all  degrees  of  freedom  except  the 
bond  variable.  This  introduces  an  extended  OQA  theory  with  one  or  several  degrees  of 
freedom  constrained  so  the  projection  of  the  solvent  modes  will  introduce  linear  couplings 
and  thus  identify  the  origin  of  the  dynamic  friction  in  the  DNM  picture.  For  the  present 
purposes,  the  formulation  will  be  confined  to  the  specific  situation  of  a  nonrotational  rigid 
bond.  From  the  perspective  of  the  OQA,  a  flexible  bond  would  allow  for  an  optimization 
of  the  bond  frequency  and  Gaussian  fluctuations  of  the  bond  variables,  whereas  the  rigid 
bond  imposes  constraints  on  the  OQA  equations.  It  should  be  noted,  however,  that  m  the 
high  frequency  limit  the  matrix  element  corresponding  to  the  bond  length  variable  becomes 
decoupled  from  other  elements  in  both  the  force  constant  and  Gaussian  width  matrices, 
so  it  actually  makes  little  difference  whether  the  bond  length  is  held  fixed  or  allowed  to 
oscillate.  Furthermore,  when  the  bond  does  not  rotate,  the  variables  corresponding  to  the 
bond  rotation  will  be  constrained  in  the  optimization.  Similarly,  when  the  center  of  the 
bond  does  not  move,  the  variables  corresponding  to  the  bond  translational  motion  will  be 
constrained. 
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as 


For  a  rigid  bond  at  rest,  the  OQA  equation  can  be  written  under  the  imposed  constraints 


{VV)g  =  0  (3-7) 

(V  :  Vy)g  =  K  ,  (3-8) 

where  the  gradients  are  in  Cartesian  space  and  the  symbol  denotes  a  Gaussian  average 
in  the  optimized  solvent  normal  modes.  In  both  cases,  the  bond  degree  of  freedom  is  appro¬ 
priately  constrained.  After  performing  the  optimization  to  find  the  ONM’s,  the  equation  for 
the  coupling  constants  is  giving  by 

(ViV,y),  =  Ci  (3-9) 

where  the  subscript  “i”  stands  for  the  ONM  mode  and  q  stands  for  the  the  bond  variable. 
Now  comes  an  important  point:  To  calculate  the  dynamical  bond  friction  in  the  DNM 
picture,  the  same  exponentially  decaying  function  f{t)  is  used  as  for  the  pure  solvent,  giving 

(()  =  .  (3.10) 

^  TT  J  U  j=l 

where  A  is  the  decay  parameter  from  the  pure  solvent  self-diflFusion  constant  [cf.  Eq.  (2.19)]. 

The  effective  DNM  spectral  density  function  for  the  friction  should  be  broadened  according 

to  the  convolution  relation  in  Eq.  (2.16),  i.e., 

IV.  APPLICATIONS 
A.  A  Simple  Example 

To  test  the  concepts  proposed  in  this  paper,  numerical  calculations  were  performed  for 
a  one-dimensional  double  well  potential,  given  by 

y{(i)  =  ^ 
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with  the  parameters  c  =  0.01,  g  =  0.l,m=  1.0,  and  p  =  5.0.  The  barrier  for  this  potential 
is  located  at  the  origin,  separating  the  two  asymmetric  wells.  Obviously,  there  are  two  sets 
of  OQA  solutions  corresponding  to  the  left-hand-side  and  the  right-hand-side  of  the  barrier. 
However,  as  the  temperature  is  increased  or  the  barrier  height  is  sufficiently  decreased, 
the  two  solutions  merge  into  a  single  quadratic  well  near  the  origin,  indicating  that  the 
thermal  excitations  overwhelm  the  barrier.  (In  contrast,  for  a  many-body  potential  there 
exists  a  complicated  branching  and  merging  of  the  multi-dimensional  OQA  solutions  as  the 
temperature  changes.)  For  comparison,  the  instantaneous  normal  mode  (INM)  correlation 
function  was  also  calculated  for  this  simple  example.  Moreover,  in  order  to  demonstrate 
the  importance  of  self-consistently  adjusting  the  equilibrium  position  to  the  center  of  the 
thermal  excitation  along  with  the  fluctuation  frequency,  the  normal  mode  equation  Eq.  (2.7) 
was  also  solved  at  the  quenched  potential  minimum  without  optimizing  the  equilibrium 
position  via  Eq.  (2.6).  The  resulting  quenched  normal  mode  spectrum  (QNM)  reflects  the 
infinitesimal  vibrations  corresponding  to  the  minima  of  Stillinger  and  Weber  s  mechanical 
stable  structures  (i.e.,  at  zero  temperature).  The  correlation  function  Cqjvm(Q  resulting 
from  the  QNM  spectrum  was  assumed  to  also  incorporate  the  exponential  damping  function. 

Monte  Carlo  importance  sampling  was  employed  to  generate  the  instantaneous  configu¬ 
rations  and  the  normal  mode  analysis  was  applied  at  each  independent  configuration;  The 
INM  frequency  was  determined  by  the  local  curvature,  while  the  ONM  equilibrium  position 
and  frequencies  were  chosen  from  the  two  sets  of  OQA  solutions  depending  on  the  instan¬ 
taneous  position,  as  was  the  QNM  frequency.  The  frequencies  were  accumulated  to  yield 
the  corresponding  normal  mode  spectra.  To  predict  the  velocity  autocorrelation  function, 
an  exponential  decay  function  with  A  =  0.1  was  employed  in  Eq.  (2.17).  For  comparison, 
the  TST  barrier  crossing  rate  for  the  double  well  was  evaluated  to  be  about  0.15,  which  is 
somewhat  larger  than  the  optimal  decay  rate.  Considering  that  the  TST  rate  in  this  simple 
well  certainly  overestimates  the  true  barrier  crossing  rate,  the  choice  of  A  is  consistent  with 
the  interpretation  that  the  incoherence  of  the  ONM’s  arises  from  the  barrier  crossings. 

In  Fig.  1,  the  velocity  autocorrelation  functions  are  plotted  for  the  exact  MD  simulation 
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results,  for  the  DNM  prediction,  for  the  QNM  prediction,  and  for  the  stable-mode  INM 
prediction.  Clearly,  the  DNM  correlation  function  gives  the  best  agreement  with  the  MD 
result,  the  QNM  correlation  function  is  out  of  phase,  while  the  INM  correlation  function 
dephases  too  quickly  after  the  second  period.  It  should  be  noted  that  in  this  example,  as 
well  as  in  the  following  two,  the  correlation  functions  calculated  by  the  different  methods 
(i.e.,  MD,  DNM,  INM,  QNM)  have  all  been  normalized  to  give  the  same  initial  {t  =  0)  value. 
Since  the  initial  value  of  a  correlation  function  can  be  calculated  exactly  from  equilibrium 
properties  through  a  Monte  Carlo  or  MD  simulation,  any  approximate  theory  can  always 
be  calibrated  to  give  the  exact  zero-time  value.  The  important  comparison  to  make  here  is 
in  the  time-dependence  of  the  correlation  functions.  It  should  also  be  noted  that  damping 
could  be  included  in  the  stable  mode  INM  correlation  function  as  in  Ref.  [10],  but  it  would 
need  to  have  a  significant  frequency  dependence  in  order  to  bring  the  INM  result  into  better 
agreement  with  the  exact  one.  That  is,  a  simple  exponential  or  Gaussian  damping  function 
in  the  INM  theory  would  just  worsen  the  agreement  between  the  exact  and  INM  correlation 
functions  since  the  stable  mode  INM  function  already  decays  much  too  rapidly  in  this  case. 

B.  Velocity  Autocorrelation  Functions  in  Liquids 

The  DNM  theory  was  next  applied  to  a  simple  homogeneous  liquid  of  particles  interacting 
through  a  pairwise  potential,  given  by 

ViTii)  =  ,  (4-2) 

where  all  quantities  such  as  mass,  length,  time,  energy  and  temperature  are  assumed  to 
be  unity.  The  numerical  studies  were  performed  at  a  temperature  of  1.2  and  a  density 
of  0.84.  After  the  system  was  relaxed  to  equilibrium,  independent  liquid  configurations 
following  every  1000  Monte  Carlo  moves  were  used  for  the  optimization.  Following  the  steps 
described  in  the  text,  liquid  configurations  were  sampled,  quenched,  and  optimized.  The 
ONM  distribution  function  was  accumulated  over  300  independent  liquid  configurations. 
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For  comparison,  MD  simulations  were  performed  for  the  same  system,  with  10^  trajectories 
being  integrated  to  yield  the  velocity  time  correlation  function. 

It  can  be  time  consuming  to  solve  the  self-consistent  OQA  equations  for  a  many-body 
system.  Fortunately,  the  thermal  fluctuation  matrix  C  is  a  relatively  small  quantity  for 
many  cases.  Such  a  narrow  Gaussian  width  allows  one  to  Taylor  expand  the  OQA  equations 
through  leading  order,  giving 

V  :  i V  :  V(V  :  V)V^C  =  K  (4.3) 

and 

vy  +  i V(V  :  V)yc  =  V(vy)  •  (q  -  q^)  (4.4) 

where  all  quantities  are  evaluated  at  the  current  optimized  position  qc.  The  above  equations 
can  be  solved  iteratively  until  convergence  is  reached. 

In  Fig.  2,  the  normalized  velocity  correlation  function  calculated  from  the  DNM  analysis 
is  plotted  along  with  the  MD  simulation  result  and  the  stable-mode  approximation  for  the 
INM  correlation  function.  [19]  An  exponentially  decay  function  with  A  =  5.0  is  used  for 
the  DNM  (I)  correlation  function,  while  a  Gaussian  decay  function  with  the  width  factor  k 
determined  from  Eq.  (2.22)  is  used  for  the  DNM  (II)  correlation  function.  It  can  be  seen 
from  Fig.  2  that,  as  expected,  both  the  INM  and  DNM  (II)  correlation  functions  agree 
with  the  MD  simulation  result  at  short  times.  At  relatively  long  times,  however,  the  INM 
correlation  function  becomes  out  of  phase  and  the  DNM  (II)  correlation  function  decays  too 
rapidly.  Overall,  it  is  thus  seen  that  the  DNM  (I)  correlation  function  with  the  exponential 
damping  function  gives  the  best  prediction  of  the  liquid  state  dynamics.  This  example 
clearly  demonstrates  the  feasibility  and  accuracy  of  the  DNM  theory  for  realistic  systems. 
In  Fig.  3,  the  exact  and  DNM  (I)  correlation  functions  are  plotted  along  with  the  ONM 
correlation  function  having  no  damping  factor  (i.e.,  A  =  0).  The  “oscillation”  is  correct 
in  the  latter  case  and  there  is  some  degree  of  dephasing  due  to  the  superposition  of  the 
different  metastable  well  solutions,  but  the  damping  function  is  obviously  required  in  order 
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to  obtain  quantitative  agreement  with  the  MD  result.  Recall  that  the  DNM  theory  explicitly 
sepcratss  the  oscillations  of  the  various  inherent  structures  from  the  damping  behavior  due 
to  the  transitions  between  such  structures. 

As  in  the  previous  example,  the  accuracy  of  the  DNM  result  illustrates  the  value  of  the 
variational  determination  of  the  ONMs.  In  particular,  the  phase  of  the  correlation  function 
oscillation  is  well  reproduced  in  the  DNM  theory  because  the  variational  effective  oscillator 
frequencies  are  chosen  to  model  the  anharmonic  thermal  fluctuations  of  the  inherent  struc¬ 
tures.  In  contrast,  the  stable  mode  INM  theory  does  not  incorporate  such  a  procedure  since 
the  stable  INMs  are  determined  from  the  instantaneous  configurations  of  the  liquid.  The 
INM  correlation  function  is  therefore  less  accurate  in  reproducing  the  phase  of  the  exact 
result,  although  a  frequency-dependent  damping  function  [10]  might  improve  the  accuracy 
of  the  INM  theory  (and,  of  course,  the  accuracy  of  the  DNM  theory  as  well).  Unfortunately, 
the  determination  of  such  a  function  from  first-principles  or  otherwise  is  not  straightforward. 

C.  Dynamical  Friction  on  Solute  Bonds 

In  this  example,  the  system  was  the  same  as  the  liquid  described  in  the  previous  subsec¬ 
tion  except  that  two  of  the  atoms  are  not  allowed  to  move.  As  was  outlined  in  Sec.  Ill,  the 
solute  molecule  was  rigid  and  held  fixed  with  a  separation  of  unit  length.  The  solvent-solvent 
and  solute-solvent  site-site  interactions  were  again  given  by  the  repulsive  1/r^^  potential. 
To  solve  the  self-consistent  OQA  equations,  the  Gaussian  average  was  expanded  through 
second-order  and  the  solution  was  found  iteratively.  The  details  are  very  similar  to  those 
described  in  Sec.  IVB.  In  the  presence  of  the  rigid  solute,  however,  the  solvent  spectrum  was 
modified  accordingly  due  to  the  presence  of  the  solute.  In  particular,  the  three  translational 
invariants  were  broken,  which  correspond  to  three  nonzero  frequency  normal  modes.  The 
dynamical  friction  on  the  bond  in  the  DNM  theory  was  then  given  by  Eq.  (3.4)  with  the 
exponential  decay  function  and  decay  constant  taken  from  the  pure  liquid  described  in  Sec. 
IVB. 
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In  the  exact  MD  calculation,  the  solvent  force  parallel  to  the  rigid  bond  direction  fu 
was  projected  out  at  each  time  step,  giving 

F{t)  =  \{Fx{t)-F2{t))-fu  (4-5) 

where  Fi{t)  is  the  force  on  atom  i  at  time  t.  [16]  The  factor  of  1/2  arises  because  the 
mass  associated  with  the  coordinate  is  the  reduced  mass  of  the  diatomic  bond.  The  force 
autocorrelation  function  was  then  averaged  over  MD  trajectories  to  give  the  dynamical 
friction  kernel  ri{t)  of  Eq.  (3.6). 

In  Fig.  4,  the  dynamical  friction  kernels  are  shown  as  calculated  from  the  DNM  analysis 
with  the  damping  parameter  taken  from  the  pure  solvent,  the  INM  stable  mode  approxima¬ 
tion,  and  the  MD  simulation.  All  curves  are  normalized  to  be  unity  at  t  =  0.  Again,  good 
agreement  between  the  DNM  and  exact  friction  kernels  is  obtained,  confirming  that  the 
DNM  model  can  be  an  accurate  approach  for  the  calculation  of  dynamic  friction  on  solute 
bonds.  The  stable  mode  INM  approach  is  again  less  accurate  in  this  exxample,  perhaps 
requiring  some  kind  of  frequency-dependent  damping  function  [10]  to  improve  its  agreement 
with  the  exact  result. 


V.  CONCLUDING  REMARKS 

In  this  paper,  a  rigorous  definition  of  the  inherent  liquid  state  structures  and  their 
metastable  normal  modes  of  vibration  was  developed  in  order  to  calculate  liquid  state  corre¬ 
lation  functions.  Prom  this  perspective,  equilibrium  and  transport  properties  can  be  studied 
in  a  systematic  fashion.  Though  the  exponential  decay  assumption  for  the  optimized  nor¬ 
mal  modes  awaits  a  rigorous  derivation,  the  intuitive  picture  of  the  damped  oscillators  is  a 
compelling  one  which  seems  consistent  with  the  linear  regression  hypothesis.  Furthermore, 
the  damping  factor  from  the  DNM  solution  for  the  pure  solvent  can  be  used  along  with  the 
OQA  theory  to  predict  the  friction  on  solute  motions,  thus  providing  a  microscopic  theory 
for  the  GLE  friction  kernel. 
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The  OQA  and  DNM  equations  can  be  applied  to  a  wide  range  of  problems.  As  an 
example,  the  activated  barrier  crossing  problem  in  condensed  phases  can  be  treated  as  an 
effective  multidimensional  quadratic  system  coupled  to  an  unstable  degree  of  freedom  so 
that  the  standard  TST  approach  can  be  used.  [13,32]  The  Gaussian  bath,  the  unstable 
mode,  and  the  linear  couplings  can  be  solved  from  the  extended  ONM  equations  with  the 
unstable  coordinate  constrained  at  the  barrier  top.  This  procedure  leads  to  a  Kramers- 
Grote-Hynes  type  rate  constant,  [32-34]  but  it  incorporates  in  a  self-consistent  fashion  the 
anharmonicity  in  the  vicinity  of  the  barrier,  the  nonlinearity  of  the  bath,  and  the  nonlinear 
couplings  between  the  bath  and  the  reactive  coordinate.  This  theory  can  also  be  extended  to 
the  quantum  mechanical  limit,  improving  upon  a  result  derived  previously  by  one  of  us.  [35] 
These  and  other  applications  of  the  present  theory  will  be  the  subject  of  future  publications. 
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APPENDIX  A:  QUANTUM  DAMPED  NORMAL  MODE  THEORY  OF  LIQUIDS 


In  the  quantum  mechanical  limit,  the  OQA  equations  are  written  as  [11] 

(VU(q  +  q))c  =  0 
(V:  VU(q  +  q))c  =  K 


(Al) 

(A2) 


where  K  is  the  optimized  effective  force  constant  matrix  and  where  V  is  the  partial  derivative 
vector  Vi  =  d/dqi.  The  notation  (•  •  •)c  here  denotes  a  multidimensional  Gaussian  average 
centered  at  q,  i.e., 

(U(q  +  q))c=  /  ^  -  /dq  V(q  +  q)  exp[-q-C~^-q/2]  .  (A3) 

Wdet[27rC]  ■’ 

The  Gaussian  width  factor  matrix  C,  in  this  case,  can  be  formally  expressed  as 

OO  1 

C  =  i:  [/3mnJ  +  /3K]“  ,  -  (A4) 


where  m  is  the  3iV— dimensional  mass  matrix  and  Qn  —  27r77./^/3.  A  unitary  matrix  U  can 
be  found  to  diagonalize  the  mass-scaled  force  constant  matrix  K,  giving  the  quantum  ONM 
frequencies 


U^KU=[I-5^]  , 


(A5) 


where  {(Di}  is  the  set  of  the  eigenfrequencies  and  [I  •  denotes  a  diagonal  matrix  with  the 
diagonal  element  given  by  u'f.  The  Gaussian  width  factor  matrix  in  Eq.  (A4)  can  be 
determined  from  the  relation 


C  =  U  [l-a]  U^ 


(A6) 


where  U  =  and  the  individual  elements  of  the  normal  mode  thermal  width  vector 

are  given  by 


1  f  {npQi/2)  ] 

PCuf  \tanh{hPQi/2)  j 


(A7) 
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Thus,  the  set  of  optimized  frequencies  {(!)}  and  average  positions  {q}  are  variationally  ob¬ 
tained  as  the  self-consistent  solution  to  the  transcendental  matrix  Eqs.  (Al)  and  (A2)  in 
iV-dimensional  space. 

One  can  next  define  a  quantum  ONM  spectrum,  giving 

1  3N 

DoNmI^)  =  ~  ^i(Qo)])qo  ) 

where  cI’i(Qo)  eigensolutions  to  Eq.  (A5)  for  the  region  R{1)  mapped  from 

an  instantaneous  liquid  configuration  Qq.  Following  the  DNM  prescription,  the  quantum 
velocity  correlation  function  for  a  simple  atomic  fluid  is  given  by 


CoNM{t)  —  J 

where  the  quantum  mode-weighting  factor  is  given  by 


(A9) 


{hpu/2) 

tanh(/i)da;/2) 


(AlO) 


It  should  be  noted  that  the  quantum  generalization  of  the  DNM  theory  is  particularly 
significant  because  information  on  quantum  dynamics  is  very  difficult  to  obtain  for  many- 
body  systems  using  direct  computer  simulation  techniques  (as  opposed  to  the  classical  case) . 
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FIGURES 

FIG.  1.  A  plot  of  the  velocity  autocorrelation  function  for  the  classical  potential  given  by 
Eq.  (4.1).  The  solid  circles  are  the  MD  simulation  results,  the  solid  line  is  the  DNM  result,  the 
dashed  line  is  the  QNM  result,  and  the  dash-dot  line  is  the  INM  result  with  the  unstable  modes 

removed. 

FIG.  2.  A  plot  of  the  velocity  autocorrelation  function  for  a  liquid  with  a  pairwise  potential 
given  by  Eq.  (4.2).  The  solid  circles  are  the  MD  simulation  results,  the  solid  line  is  the  DNM  (I) 
result  for  an  exponential  damping,  the  dashed  line  is  the  DNM  (II)  result  for  a  Gaussian  damping, 
and  the  dash-dot  line  is  the  INM  result  with  the  unstable  modes  removed. 

FIG.  3.  A  plot  of  the  velocity  autocorrelation  function  for  a  liquid  with  a  pairwise  potential 
given  by  Eq.  (4.2).  The  solid  circles  are  the  MD  simulation  results,  the  solid  line  is  the  DNM  (I) 
result  for  an  exponential  damping,  the  dashed  line  is  the  ONM  result  for  no  damping  (A  =  0),  and 
the  dash-dot  line  is  the  INM  result  with  the  unstable  modes  removed. 

FIG.  4.  A  plot  of  the  friction  kernel  for  a  rigid  solute  bond  as  described  in  Sec.  IVC.  The 
solid  line  is  the  MD  result,  the  dashed  line  is  the  DNM  prediction  with  the  exponential  damping 
parameter  determined  from  the  pure  solvent,  and  the  dash-dot  line  is  the  INM  result  with  the 
unstable  modes  removed. 
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